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Chapter 1 

Abstract 



Abstract 

Several experimental results show that it is possible to extract useful phase information 
from reflected GPS signals over the oceans. In this work we begin the development of the 
theoretical background to account for these results and fully understand the phenomena 
involved. This information will then be used to define and carry out new experiments 
to evaluate the feasibility of using the phase from reflected GPS signals for altimetric 
purposes and the advantages of using interferometric combinations of the signals at different 
frequencies — the PIP concept. 

We focus on the coherence properties of the signals, including the PIP interferometric 
combination of phases in the different frequencies. In this work we will concentrate on a 
static, 8 m high receiver (at least in regards to the simulations), and an infinitely removed 
static source. As the ocean moves, the received field will pick up a random phase. We 
want to understand the behavior of this phase, as the goal is to carry out altimetric mea- 
surements using phase ranging. We will also keep in mind that this random phase carries 
geophysical information (intuitively, the bigger the significant wave height, the larger the 
phase excursions). 

Our simulations are based on the Fresnel integral and use simulated Gaussian oceans 
using the Elfouhaily et al. spectrum. The simulation tool, FRESNEL, is capable of gen- 
erating time series of the reflected field at different frequencies, and then analyzing their 
properties. This software is written in IDL. 

The most important point we need to answer is whether the signal can be tracked in the 
aforementioned situations. As we show, the PIP combination of the signals helps clean the 
signal from noise — the more correlated the signals at different frequencies the more effective 
is the PIP mechanism. The following questions are specifically addressed: 

• What is the reflected field spectrum? This is determined by the orbital ocean motion, 
for a static receiver. We also discuss the moving case and show that it depends on 
the system gain. 

• Does the winding number accumulate, or does it average to zero? Winding number 
refers to phase accumulation. Large phase excursions (many cycles) are seen in our 
simulations. Our simulations show, however, that there isn't a preferred direction for 
phase winding. This means that altimetric phase measurements can be accurate. 

• What are the theoretical values for the average values of the interferometric fields? 
How do the fields correlate across frequencies? We show that the most important 
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factor is the relation between significant wave height and the (real or synthetic) elec- 
tromagnetic wavelength. 

• What is the coherence time of the signals (the reflected phase)? We show it is longest 
for the PIP combination. For rough seas, the correlation between the fields (and 
therefore interferometric coherence) in different frequencies disappears, and the co- 
herence time goes to zero even for the interferometric combinations. In calmer ocean 
conditions, however, our results indicate that the interferometric combination remains 
coherent while the individual signals lose coherence rapidly. 

• What is the structure function of the reflected phase? We see a good fit with a random 
walk model, with a drift rate proportional to wind speed. 

• In what ways is the PIP interferometric signal different and, presumably, superior to 
the original ones? Coherence, a basic element for altimetric purposes. 

Finally, we discuss the robustness of altimetric phase measurements after low-pass fil- 
tering the PIP combination. We show that low-pass filtered PIP data can provide more 
robust and accurate altimetric measurements to detect slow-varying geophysical signals in 
the ocean. 
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Chapter 2 

Introduction 



This is the sum of Technical Notes PIPAER-IEEC-TN-1100 and PIPAER-IEEC-TN-2100, 
Theoretical Feasibility Analysis. The inputs to this WP are 

1. GMV/IEEC Proposal GMVSA1123/99 (GMV) 

2. Relevant literature 
The outputs are 

1. This technical note 

2. Recommendations for experiments and post-processing procedures 
And the Tasks: 

1. Review the relevant documentation 

2. Establish the worst and best case scenarios for the concept applicability 

3. Issue recommendations for experiments and post-processing procedures 

4. Prepare this technical note 

In addition we analyze some issues relating to the architecture of the proposed PIP 
instrument. 
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Chapter 3 

State-of-the-Art, Review 



3.1 Historical overview 



The PIP idea can be traced back to the PARIS concept \ Martin- Neira, 1993\ , where the 



use of GNSS signals as sources of opportunity for bistatic altimetry was first proposed. 
The goal of the present research is to assess the possibilities of using the phase information 
of the reflected signal from a PARIS GNSS system for altimetric purposes. In particular, 
we will study the merits of using the phase from a dual frequency GNSS system in an 
interferometric fashion. Other uses of reflected phase data can also be envisioned (surely 
the temporal behavior of phase fluctuations and other characteristics contain geophysical 
information about the sea surface) but we will not explicitly investigate this possibility in 
what follows. Rather, we will, for now, content ourselves with providing the background 
for discussion of these ideas in future work. This is a natural byproduct of our interest in 
phase altimetry, as we will see. 

That the GPS signal reflected from the ocean surface can retain coherence under some 
conditions has been demonstrated in the past. In \Auber et al, 1994\ we have the first 



report of a GPS receiver "locking" onto a sea-reflected signal. Motivated by this somewhat 
surprising event, tracking of reflected signals was later seeked and achieved for extended 
periods of time in a series of experiments with flights (over the Chesapeake Bay and the 



Eastern Shore of Virginia) at an altitude of up to 5500 m Garrison et al, 1996\. This 
was accomplished using an Ashtech Z-12 off-the-shelf receiver with a LHCP nadir looking 
antenna, and it was observed that carrier lock was obtained. This particular receiver would 
not have been able to track otherwise. 

On the other hand, several groups have to date successfully analyzed reflected data 
in order to correlate the direct and reflected signals with clean code replicas. In the ex- 
periment carried out by ESA in September 1997 two receivers were employed, operating 
independently of each other, one tracking the direct signal with an up-looking antenna, the 
other recovering the reflected signal via a down-looking LHCP 9 dB helical antenna over an 
18 m bridge near Rotterdam (The Netherlands), the "Zeeland Brug". After amplification, 
each of the signals was sent to a different GPS receiver (GEC Plessey builder kit 2). The 
receivers down-converted the signals to IF (4.309 MHz, with a bandwidth of 1.9 MHz). 
The IF signals were sent directly to a high-performance sampling card, which sampled the 
data at 6.25 MHz with 2 MHz bandwidth. They were then cross-correlated with replicas of 
the signal adjusted for Doppler due to satellite motion and with 20 ms duration (20 C/A 
code periods). Using the correlation delay data from the direct and reflected signal, this 
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group solved for the height of the bridge over mean surface (which depended on the tide) 
as well as a hardware bias constant via a least squares fitting procedure. In addition, the 
receivers recorded standard RINEX files, including SNR. The GPS signals were deemed 
usable if the SNR of the reflected signal was greater than 6 to 9 dB, and if the geome- 
try was favorable (visibility and multi-path were rather complicated in some geometrical 
configurations, which were then excluded). Significant Wave Height during the experiment 
was around 1.4 meters (equivalent to a height standard deviation of about 35 cm), and 
the altimetric performance was rather poor, as only the C/A code could be used for corre- 
lation (in general, a 1% of the chip length is assumed to be the optimal precision for one 
measurement — this is 3 m for C/A code). A very good description and further details can 
be found in \ Caparrini, 199^ . 

The IEEC Earth Sciences Department has been able to reproduce the analysis of the 
Bridge Experiment data with similar results. The analysis tools for the generation of the 
GPS signals and their correlation with the received ones have been developed in Matlab. 

Another aspect of the Bridge Experiment analysis carried out by ESA was the attempt 
to extract phase information from the data. The first attempt of phase processing was an 
evaluation of Doppler frequency based on zero-crossings counting \ Caparrini, 199^ . The 
approach was to first multiply the signal by an appropriate clear replica of the PRN code 
of a satellite in view. It was then expected that a peak in the Fourier transform of the 
resulting signal would appear at the IF plus the Doppler frequency. The direct approach 
was not deemed feasible, however, due to the limited length of the signals (only 10 ms long), 
which did not allow for enough spectral resolution. The signal multiplied by the replica was 
instead first filtered with a rather narrow filter in order to eliminate all but the expected 
Doppler-shifted component. After that, frequency was estimated using the zero-crossings 
count for each 10 ms stream of data, with rather good results for both the direct and the 
reflected signal. This suggested that the reflected signal still contained a certain degree of 
coherence. A second study was then devoted to the measurement of the relative tide height 
(that is, the difference in the height measurements at different times through the difference 
in the number of zero-crossings between the direct and the reflected signals). The reasoning 
was that if the phase could be measured at all, it would be only up to a phase bias. If the 
bias is assumed to be constant over some short measurement times, then the change in the 
phase between those times would be related directly to the change in delay between the 
direct and the reflected signals. This calculation, however, could not be carried through 
because the data was too noisy and the intervals between data measurements too far apart 
in time. 



3.2 The goal of this research 

As we have seen, several experimental results show that it is possible to extract useful phase 
information from reflected signals. In the analysis sketched below we will begin to develop 
the theoretical background to account for these results. The goal of this work is to continue 
this development and analysis of past work to fully understand the phenomena involved. 
This information will then be used to define and carry out new experiments to evaluate the 
feasibility of using the phase from reflected GPS signals for altimetric purposes. 

We will focus on the coherence of the signals, including the interferometric combination 
of phases in the different frequencies. In this work we will concentrate on a static, 8 m high 
receiver (at least in regards to the simulations), and an infinitely removed static source. As 
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the ocean moves, the received field will pick up a random phase. We want to understand 
the behavior of this phase, as the goal is to carry out altimetric measurements using phase 
ranging. We will also keep in mind that this random phase carries geophysical information 
(intuitively, the bigger the significant wave height, the larger the phase excursions). 

The PIP concept concept consists in combining two phases from the signals at different 
frequencies. Imagine a very smooth and slowly changing ocean. By combining the two 
phases we obtain a slowly changing interferometric phase. In some sense we would expect 
that this interferometric phase will be less sensitive to small ripples in the ocean and thus 
provide a more robust ranging tool for altimetry. 

Several implementations of this concept can be envisioned: a static receiver on a coast or 
bridge to monitor tides or floods as well as sea state. Or deployment on boats, aircraft and 
spacecraft for global altimetric or oceanographic measurements. Carrier altimetric signals 
over the oceans, even when seen by a fast moving Low Earth Orbiter, will mostly be low 
frequency. A strong topographic signal over the oceans would be a slope of 1 meter every 
100 km. As seen by the LEO (traveling at 7 km/s, say), this translates into 0.35 Hz in LI, 
a rise of 5 cycles in 100/7 seconds. In L25, this becomes 0.01 Hz. 

On the other hand, the roughness and motion of the sea will induce high frequency 
"jitter" on the received signal. This is a nuisance for altimetry, and this is where we expect 
the PIP concept to bring added value, as the interferometric combination should be less 
sensitive to these effects. It is appropriate to recall here the old ado, 'What is noise to 
some is signal to others." This jitter may contain very useful information. 

In the author's mind, the most important point we need to answer is whether the signal 
can be tracked in the aforementioned situations. As we will see, the PIP concept will help 
clean the signal from noise. For the most part, we will focus in this work on the low altitude, 
static receiver situation. The moving cases at higher altitudes will be left for future work, 
as we will discuss in the last Chapter. 

In this report we will address, among others, the following questions: 

1. What is the reflected field spectrum? 

2. Does the winding number (to be described below) accumulate, or does it average to 
zero? 

3. What are the theoretical values for the average values of the interferometric fields? 
That is, how do the fields correlate across frequencies? 

4. What is the coherence time of the signals (the reflected phase)? 

5. What is the structure function of the reflected phase? 

6. In what ways is the interferometric signal different and, presumably, superior to the 
original ones? 
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Chapter 4 



Background and Applicable 
Models 



The Rayleigh criterion defines a surface to be smooth if < X/ (8 cos 9), where is the 
surface height standard deviation, 9 is the incidence angle with respect to the normal to 
the surface. A more stringent condition is provided by the Fraunhofer criterion, which is 
used to define the far-field distance of an antenna, < A/(32cos0). A surface becomes 
smooth under two conditions: cr^ ~ or 9 ~ 90°. The effective roughness of the surface is 
therefore cjq cos 9 \ Beckmann et al., 196% . 



Several aspects must be considered to understand the scattering phenomenon. We will 
examine two limiting situation models. In one case, we can assume that the surface is 
made up of many independent mirrors (specular points) whose contributions to the field 
add up incoherently. The other limit is that of a surface which deviates from flatness only 
slightly. These two cases must be treated differently. The goal in both cases will be to 
understand the behavior of the received phase and modulus. Let us keep in mind in the 
back of our heads that, at the end, we should look at the relative phase between the two 
available wavelengths, which should be a statistically better behaved quantity than either 
of the two phases. 

4.1 Statistics of the field 



In this subsection we briefly review some of the results in [ Beckmann et al, 1963^ , concern- 
ing the sum 



n 



U = re # = £y^- (4.1) 

3=1 

This is relevant to understand what happens when we sum the fields in a rough reflection 
situation. We state the result for two illuminating conditions: 

1. A uniformly distributed phase from —it to ir. 

2. A normal phase distribution with standard deviation a. 

In the first case the resulting phase has again a uniformly distributed phase, but the modulus 
has the so-called Rayleigh distribution p(r) = ^ exp(— r 2 /n) — and therefore has a non-zero 
average (the average is proportional to n 1 / 2 ). In the second case we find that the components 
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of the resulting phasor are normally distributed (which is a quite general result following 
directly from the Central Limit Theorem) and that 



(r 2 ) = n 2 e- a + n(l - e~ a ). (4.2) 

In this beautiful expression we can see the coherent and incoherent contributions to the 
field modulus. Beckmann summarizes: 

Outside a narrow cone (or wedge) about the direction of specular reflection, 
the amplitude of the field scattered by a rough surface is always Rayleigh- 
distributed; if the surface is very rough, and grazing incidence is excluded, the 
amplitude of the scattered field is Rayleigh-distributed everywhere. 

The specularly scattered field is composed of a coherent component and a random, Hoyt- 
distributed component. When the surface is very rough, the latter becomes incoherent and 
the former vanishes. In fact, if the surface height distribution is normal with deviation a^, 
then 

( r 2 ) = n 2 e -(4™ c co^/A) 2 +n(1 _ e -(^a c cose/Xfy ^ 

Thus, we see that the magnitude of the reflected field should depend mainly on the ratio 
of significant wave height to wavelength. We will obtain a related result through more 
sophisticated analysis below. 



Some verification of these ideas via simulations can be found in \ Daout et al, 199i ], 
where the authors, using a 1-D bistatic scattering simulator, reproduce the mentioned 
aspects of the phase statistical distribution. The surface is there simulated by facets, and 
the Kirchhoff theory for the scattering field from each facet is used to accumulate the rays 
with regard to their phase. The phase histograms they show verify very clearly the ideas 
just discussed. In particular, that the signal is very coherent at low elevations. 

4.2 The Fresnel-Huygens-Kirchhoff integral for the field 

This is basically the scalar Kirchhoff approximation, valid for surfaces with large radii of 
curvature compared to wavelength. Let the incoming field be described byQ 

pikr 

Uo(p) = • (4.4) 



The Fresnel integral for the scattered field is (see \ Born & Wolf, 199% p. 380, eq. 17]) 



_„■ f „ik(r+s) 

U(p) = — / n (q- n) dS. (4.5) 

47r J rs 

The vector q = (q±, q z ) is the scattering vector, the vector normal to the plane that would 
specularly reflect the rays in the direction we are looking. This vector is a function of the 
incoming and outgoing unit vectors Hi and n s , q = k(hi — h s ). The scattering vector is 
related to the specular angle through cos f3 = z ■ q/q. Note, as an aside, that in the 
nadir case all wind direction dependence disappears (surface anisotropy). Changing wind 
direction is akin to performing a rotation in the surface, and the integral is invariant under 
rotations in the nadir case. The nadir case is special because R[q(x)] = q(R[x]). 

Here we do a scalar treatment, think of U as a component of the field. 
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4.3 Geometrical Optics 



This is the process that dominates in specular scattering when the frequency is large (i.e., 
the wavelength is small compared to the wave height and to the surface correlation length) . 
The case of L-band scattering from the ocean is at best borderline "high-frequency" , and 
we are carrying out some separate studies within the framework of GNSS-OPPSCATQ to 
understand the relevance of this approximation in such circumstances. To this end, we have 
developed the speckles software tool, written in IDL (Interactive Data Language). In Ge- 
ometrical Optics the surface height standard deviation is assumed to be at least of the order 
of the electromagnetic wavelength. According to the analysis in [Beckmann et al, 1965], 



the resulting wave will be largely incoherent. There is an important ingredient in our sit- 
uation, however, and that is because we have a Woodward Ambiguity Function (WAF) at 
our disposal to select a given surface patch, as we will see shortly. This means that we can 
filter field contributions from a rather small area. 

In the Geometrical Optics approximation to the Kirchhoff theory for electromagnetic 
scattering, the physical picture can be understood in terms of a specular point model. 
That is, the field at the receiver is the superposition of the fields generated by a number of 
"mirrors" (not flat mirrors, thought, rather parabolic caps) on the scattering surface which 
are oriented in the correct manner. The radiation from each specular point is as coherent 
as the incoming radiation, but there is no coherence in the phase relationship between 
the radiation out-coming from the different mirrors if the surface height distribution has 
a large range. The result is that there is power at the receiver, but that the phase of the 
received signal is randomly distributed and the power is simply proportional to the number 
of scatterers. The behavior of the resulting phase will depend on the geometry, on the 
number and distribution of mirrors, and also on the temporal variation of these quantities. 
The goal of this section is to understand this relationship and quantify it. 

It is useful here to stop and give a physical picture of the specular point model — one of 
the two ingredients in Geometrical Optics. We will treat here the radiation coming from 
each mirror separately, as was mentioned above. Since the size of each mirror is small 
compared to the Fresnel zone, we can use the Fraunhofer approximation (plane waves all 
the way through) to estimate the out-coming field. 

Recall that the Fresnel integral for the scattered field is 

—j (■ e ik(r+s) 

U{p) = — Tl (q- n) dS. (4.6) 

4-7T J rs 

In the far field (both emitter and receiver are very far compared to the size of the scatterer, 
i.e, the size of the scatterer is much smaller that the Fresnel zone), and since our integral 
is near a specular point with n ~ q, we find that 



-i e 



ik(r'+s') 



U(p) = — / K ■ e~^ r (q- n) dS, (4.7) 

47r r s 



see Figure 4.1 for a pictorial definition of the involved vectors. Here all we have done is 
to approximate k{r + s) « k(r' + s') — q ■ r. Now, let us focus on a single specular point: 
imagine that there is only one specular point on the surface. This means q ■ n = q. To 
perform the integral let us use a coordinate system in which the z axis is parallel to q. In 



2 ESA Contract 13461/99/NL/GD, Utilization of Scatterometry Using Sources of Opportunity. 
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z 




Figure 4.1: Vector geometry for scattering. 



this coordinate system the tangent plane at the specular plane is therefore parallel to the 
x-y plane. Hence, dS = 4-f = d 2 x. We need to compute 



I=qjn-e~ iq < x ' y) d 2 x, (4. 



where 



■ ik(r'+s') 

u = ir^^ L ( 4 - 9 ) 

47r r's' 

Now we use the stationary phase approximation. The idea is that as q gets larger the 

integral gets contributions only very near the specular point — the contributions farther out 
cancel out. Then, 

/ « I = qTZ spec e- iqZs ^ J e"^ ^ z ^ XiXj (fx (4.10) 

2-7T 

= qK spec e-*« r ^ r (4.11) 

-iq det 1 / 2 ( V?-z 



spec 



i 7^- S p ec 6 $ p . (4.12) 

det 1 / 2 (Viz ) 

\ J spec J 

This result, that the square root of an inverse determinant follows from a Gaussian integral 
is a much used result in Quantum Field Theory (see for instance, Quantum Field Theory 
by L. H Ryder, Cambridge Univ Press 1996 for a derivation). It can be seen to follow from 
extending 

to higher dimensions. Now, it is not hard to show that det 1 / 2 ( V 2 z ) = Jr\T2 — 

\ J spec J 

the determinant just yields the products of the radii of curvature at the specular point: 
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the determinant is invariant under coordinate transformations of the surface, so we can 
use a coordinate system along the principal directions — the surface is then described as 
a simple parabola. The result is then immediate (the radius of curvature is the inverse 
of the second derivative in such coordinates and the matrix becomes diagonal in such a 
coordinate system, so the computation of the determinant just yields the product of the 
diagonal terms). Finally, 



-i e 



ik{ r'+s') 2TT 



U (p) = - A ry— ■ in spec e- tq - r ^ (4.13) 



47r r ' s ' ' detWv? 



ij z 



spec 



& ^KlVspec~\~S spec ) 

-y/rjfi (4.14) 

The coefficient 7Z spec depends on the Fresnel coefficients and on the local geometry — we 
will discuss this below. The cross section is given by the ratio of the resulting field squared 
divided by the field squared at the surface (given above) times Anr 2 . The result is 



a = ' 2 ' = \H spec \ 2 7r ri r2. (4.15) 



For all this to work, some requirements have to be met on top of those from the Kirchhoff 
approximation. Basically, this is a high frequency approximation — q has to be large. Let 
us look at things in ID, along one of the radii of curvature. The integral is essentially of 
the form 

J = [** dxe iqx2 ^ (4.16) 
J-u 



(4.17) 




For the integral to become ^/tt, we need l x ^q/r x to be large — of the order of 10. The 
expression for the incoherent sum of all the specular points is now immediate from Equation 



4.14, given a large height deviation. This expression for the field is also useful for a not-so- 



incoherent sum of the fields — we will return to this point below. 

4.3.1 Corrections to the Geometrical Optics approximation: frequency 
dependence 

In this section we add a comment on the next order corrections to the Geometrical Optics 
approximation. This is not a crucial section for the work at hand and can be taken as a 
small aside, but we believe it will be useful in future work. 

In reality the surface is not a parabola and other terms appear. These are frequency 
dependent corrections to the field. In order to obtain the corrections we need to expand 
the integrand as a Taylor series beyond the well-behaved parabolic term. The basic idea 
for this type of computation is that 



1(b) = j e- ax2+bx dx = e~ b2 ' Aa ^Ta, 

and that 

in 

e- ax2+hx x n dx. 



d n 



b=0 
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From this we can compute 
The expansion is 



e- ax2 x 4 dx 12 



I = qJ-R- e- iqz(x ' y) d 2 x = I + I 4 + . 



where 



Now, 



Ta — aV p- ic i z spec lq V 4 ., ,z 

14 — q l\. spec ti ^ V yfc/^ 



spec 



x l x J x'"x"e 



1. 1 —i% V 2 -z\ XiX-j ,9 



Q i*-spec e ~ v^l-iZ 



4! 



ijkl A 



spec 



spec*- 



I'Q %spec _ 



spec 



-1 



+ 



spec 



spec 



2ir 



-iq 



det 1 / 2 ( 



V 2 -z 
13 



and 



n 



spec 1 - 



-tq Zspec _ 



q4\ 



'Aspect v jjfc;^ 



spec 



V 2 -z 

V 



spec 



-1 



spec 



(4.18) 



(4.19) 



xVxVe ia V -^Uc x ^d 2 £ (4.20) 



(4.21) 



-1 



spec j 



(4.22) 

In a future implementation of speckles (see discussion below), we would like to include 
these second order effects. These effects can also be accounted for theoretically using, 
e.g., Gaussian statistics. At the end the final expression must depend only on the usual 
parameters. Note that the I 4 has units of length. Thus, aside from constants, these 
higher order corrections may be of the form j^, say. In this last expression we see that 
this correction is indeed inversely proportional to frequency, as it should, since Geometric 
Optics is a high frequency limit. 



4.4 The WAF zone 

If the signal is filtered by cross-correlation with several C/A (or P) code periods at the 
delay and frequencies corresponding to the specular zone, for instance, only the surface 
patch in the first Delay and Doppler zone will contribute (let us call it the "WAF zone" 
to associate as well as distinguish this concept from the concept of Fresnel zone). This is 
described by the equation 

SNR(r, /.) - ^L- • ^ / x V, t, Sr, Sf) iA. 

which is a slight refinement of a similar equation in [ Zavorotny et al., 199i\ and which can 



be found in the review \ Ruffini et al, 199£ \. The important thing to keep in mind is that 



the support in the integrand, the WAF zone, is the intersection of four spatial zones: 
1. The receiver antenna footprint. 
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2. The annulus zone defined by the A 2 function in \ 2 . 

3. The Doppler zone defined by the |5| 2 function in x 2 - 

4. The scattering cross section coefficient gq. 

The WAF, through its support, selects a given portion of the surface from which power 
is measured. 

To get a feeling for these numbers, the area involved in the reflection is of the order of 
10 km 2 for an aircraft and 1000 km 2 for a Low Earth Orbiter. For h » ct, where r is 
the chip length (1 /xs for C/A code) and h the receiver height), the one- WAF (or one-chip 
here) area can be approximated by Ay/AF = Zither/ cos 2 6. Note that Awaf increases only 
linearly with altitude. 

We can add another element at this point. In the above discussion, the WAF is assumed 
to result from choosing a given delay and Doppler in the reflected signal. However, as time 
passes the signal reflected by a given surface element cannot be characterized by a fixed 
Delay and Doppler: both will change as the situation evolves. It is clear, however, that 
a more sophisticated WAF algorithm can be devised that focuses on that specific surface 
patch. This observation opens the door to longer integration times, which can be useful for 
altimetric purposes. 



C/A and P code WAF size 




100 200 300 400 500 

Height, km 

Figure 4.2: Nadir case first chip radius size for C/A and P codes. This is for a static situ- 
ation, no Doppler filtering affects the result (based on the simple formula s = \/2 * h* r c ), 
where r c is the chip length in meters (300 m for C/A code). 
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Chapter 5 



Simulation Tools: fresnel and 

SPECKLES 



These are two modeling tools we have developed to understand, via a computer simulation, 
some of the characteristics of the reflected field, fresnel models the reflected field (just the 
carrier part after removing the source time-dependence) using an ocean model and direct 
integration of the Fresnel-Huygens-Kirchhoff integral described above. We have carried out 
several types of simulations: a slowly vertically moving receiver, a static receiver, a realistic 
ocean and a so-called chaotic ocean. We will focus here on the static realistic ocean case. 

speckles has been developed with two goals in mind. We will eventually attempt 
to compute the reflected field using a GO method: finding first the specular points and 
integrating just around these to sum up the resulting field. We will just outline here 
some of its characteristics and initial results. We have found that this is very difficult to 
do (in general, comparisons with the Fresnel integral are poor). We have seen, however 
that the number of specular points correlates reasonably well with the reflected Fresnel 
field, specially in rough ocean simulated conditions. We are attempting to understand 
this phenomenon better, and to study its implications: GO, in some form of another, is 
presently the leading model used for the analysis of reflected power in GPSR. It is a simple 
method, and experimental results to date seem to validate it. We would like to understand 
its somewhat puzzling high performance, since L-band is not really high-frequency in the 
ocean surface case. That is, GO is not, strictly speaking applicable in the whole range 
of ocean conditions with an L-band instrument — the wavelength in neither large of short 
enough. Nonetheless, some semi-empirical modifications of this model are currently being 
successfully used. We would like to understand this success and, more importantly, its 
limitations. We will briefly mention some of this work, although it is not central to the 
study. 



5.1 FRESNEL 

This routine ( written in Interactive Data Language — IDL) integrates the field from the 
ocean surface below, using nadir incidence. The receiver is hovering a few meters over 
the water. The ocean model is Gaussian (we use a routine from B. Chapron using the 
Elfouhaily spectrum) \Elfouhaily et al., 1991 \, "chaotic" (random heights, moving) and a 
non-Gaussian simulations are planned for future work. The software simulates the reflected 
complex field at LI, L2 and L5 frequencies. The goals of this part of the work were to show 
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that the simulation is possible, to compare the behavior of the field in the different cases, 
and to assess the potential of the PIP idea. 



5.1.1 Random ocean 

The random ocean is generated by choosing a height standard deviation. The code assigns 
a random height with the desired characteristic to each 10 cm by 10 cm square of the 
reflecting surface. In addition, a velocity random field is chosen, and each square moves up 
and down with a chosen period. 



5.1.2 Elfouhaily spectrum 

The procedure to compute the ocean state, based on the IDL routine kindly supplied by B. 
Chapron, is, starting from the energy spectrum, to generate a plausible random spectrum. 
To do this, the square root of the energy spectrum is taken, and a uniformly random phase 
appended to it. The appropriate time dependence is also added to this random phase — 
using the dispersion relation for gravity ocean waves. Then, the inverse Fourier transform 
is taken. Using this procedure, a moving, random ocean with the appropriate spectral 
characteristics is generated. Using a uniform random phase distribution in the different 
frequencies means that the resulting ocean is Gaussian. This can be changed, of course, 
but we have not attempted to analyze the subtleties of non-uniform phase distributions in 
this work. 

Generally the sea state is characterized by its power spectrum. As the sea evolves in a 
random process, the spectral density of sea elevations is obtained by the Fourier transform 
of the autocorrelation function of the elevations C,{x,y). The spectral density can also be 
estimated by the Fourier transform of one realization of y) [Blackman and Tukey, 1958]: 

F(k x ,k y ) = \TF[((x,y)}\ 2 

If we determine the matrix <f> of random phases uniformly distributed between and 2ir, 
the sea height at the point r = (x, y) is: 



C(r)=TF- 1 [^F(k x ,k y )e^} 

Thus the probability density functions of heights and slopes are Gaussian. For this study 
we use the unified spectrum of [Elfouhaily et al, 1991]. This is just a consequence of the 



Central Limit Theorem applied to sums of harmonics. 
5.1.3 Field calculation 

The field is calculated as a complex number using the Fresnel integral expression: 

—j r „ik(r+s) 

U{p) = — \ K (q ■ n) dS. (5.1) 

47r J rs 

To be precise, we Fraunhofer-expand the part of the integrand corresponding to s the 
distance to the infinitely far transmitter. That is, we write s = ks' — k{ n ■ x in the phase 
part of the integrand, where x denotes the position of the scattering point. Thus, the 
integrand becomes, 

U(p) = / n (q- h) dS. (5.2) 



18 



We actually throw out all the s' dependence (in effect assuming an incoming plane wave, 
as well as on the Presnel coefficient ) and obtain 

U(p) = — J Tl (q ■ h) dS. (5.3) 

This means that the outgoing field would be a plane wave if the surface were flat, with unit 
modulus. 

The ocean surface is divided into squares 10 cm wide (or less), and an area of up to 
200x200 meters is integrated. For the most part we will concentrate on the first chip zone, 
with radius r = y/2 * h * r c . 

The resulting field is a complex number in this scalar treatment, U(p(t)) = re 1 ^. In 
order to evaluate the phase as the receiver moves up we take 

8<t>i = ]n(HU)/w(U-i)) (5-4) 

with w(t) = e l ^ l \ Summing the 5<fii yields the overall phase. This method is related to 
the definition of winding number of a curve around the origin in the complex plane: 

/= cf — = cf d(lnz)=J2^—- (5-5) 

Jy Z Jy . Zi—1 

A full cycle, 360 degrees, is equivalent to 19 (LI), 24 (L2), 25 (L5), 86 (L12) or 568 (L25) 
cm. We will return to this definition shortly. 



5.2 SPECKLES 



As mentioned above, this code has been developed to study the use of Geometric Optics for 
the study of ocean surface reflections in L-band. It will also be used to study specular point 
statistics, using both Gaussian and, eventually, non-Gaussian models. Taking a Gaussian 
ocean model, we determine the specular points distribution and compute the reflected field 
using geometrical optics. We are to compare this to straight integration of Fresnel integral. 
Then the analysis will be compared with a non-linear ocean model. 



5.2.1 Specular point determination 

We consider the incident radiation as a plane wave normal to a Gaussian surface and a 
receiver at (x r , y r , z r ). The point (x, y, () of the surface is considered to be a specular point 
if the slopes follow rather clear geometrical conditions (see for instance [D.E.Freund 1997]): 

dC - X ~ Xr (5.6) 



dx z r — C, 

d( y-y T 



(5.7) 



dy z r -( 

Our approach consists in the detection of slope variations around the expected specular 
slopes defined above. We compute at each point P of the surface, and for both directions 
x and y, the expected specular slope S — i.e. the slope needed to reflect the incident wave 
to the receiver — and the slopes of the facets located just before (slope S [before]) and after 
(slope S [after]) the point P. The condition for P to be a speckle is that (S[before]-S) and 
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(S -S [after]) must be of the same sign. This determination is attractive because we doesn't 
need to define any tolerance on the specular slopes. We are neglecting saddle point by this 
approch (saddle points are not well defined in discrete spaces). 

Then an analysis is made on the radius of curvature at the specular point. The Kirchhoff 
approximation leads to a condition on the radius of curvature at the specular facet, 

yjqy/wi > 10- (5.8) 
The product of the principal radii of curvature can be computed using the derivatives of 



1110= 8 m/s 




Contours at 0, 4, 10 and 20 lambdas, L1 and L2 

Figure 5.1: Fresnel contours for the 10 by 10 surface patch under and 8 meter high 
receiver. U10=8 m/s. The colored contours represent the ocean instantaneous topography, 
with higher areas represented brighter. The line pairs are iso-lambda delay curves for LI 
and L2 (4, 10 and 20 lambdas are shown). 
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C(r), see [Barrick, 1968]: 



\nr 2 \ 



(i + C 2 + C 
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yy 



C 2 

Saf 



■iy 



The following Figures ( |5.2j to 5.4 ) present the specular point position (circles) for three 
different times on the moving surface. The receiver is located at 8 meters above the surface, 
which has a 5 m by 5 m dimension and 1 cm of resolution. The stars correspond to the 
specular points that are under the curvature condition. The specular distribution is greatly 
modified within 10 ms. 




Figure 5.2: Specular points position for Time = s 




Figure 5.3: Specular points position for Time = 10 ms 



5.2.2 Field computation 

The idea is to compute the Fresnel field with the contribution of scatterers only. We 
integrate on the vicinity of each mirror. The result is then compared to the straight Fresnel 
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Figure 5.4: Specular points position for Time = 20 ms 



computation. To date these comparisons have not been very successful. It appears that 
in order for the Geometric Optics approach to work at this level (not after substantial 
averaging) we need to impose more stringent frequency or surface scale characteristics 
and/or take into account higher order effects. If the main contribution is from specular 
points, as we suspect, the key may be to calculate the field from each scatterer to higher 
order than the stationary phase approximation. Work in this area is planned for the future. 
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Chapter 6 

Coherence and Structure Functions 



6.1 The coherence and structure functions, and the coher- 
ence time 

Another question that needs to be answered is how the sea surface motion and varying 
geometry affect the coherence of the reflected signal. There are several models for the 
sea spectrum, and these can be used to try to extract this information. Consider the 
toy model of one-dimensional scattering, with static, lined up, receiver and transmitter, 
with the scatterer (a single facet) also in line but now moving up and down with average 
speed v. It is easy to see that the correlation of the direct and the scattered signal is zero 
unless the coherent integration time is less than the characteristic time of the surface at 
GPS wavelengths. The coherent integration time is roughly that time for which the RMS 



phase error is 1 radian | Thomson et al, 199C \. More precisely, let the coherence function 
be defined by 



c(T) 



T e^dt 



o 



(6.1) 



where (j){t) is the error phase. Then, the coherence time is the time T it takes for (C(T) 2 } to 
drop to 0.5 — say { Thomson et al., 199(\ . Note that for long times the coherence function 



can also be related to the magnitude of the zero-frequency spectral component of the 
(windowed) normalized field associated to the error phase (that is, in our simulations, the 
Fourier transform of the normalized field). We will get back to this point below. 

The coherence time of the GPS signal is about 2-3 hours. What about the reflected 
signal? If the WAF zone and the sea surface are fixed, we should expect the same. If we 
include time dependence (and we ignore all motions), we can rewrite the equation on the 
statistics of the signal by 

U(t) = re** = e^~ iujt = e~ iujt ^ e*** , (6.2) 

3 3 

so the resulting field is as coherent as the incident one. But in real situations neither of this 
is the case. Both the ocean surface and the WAF patch are moving (although this second 
effect can be corrected in principle by the use of S AR techniques) . The characteristic time 
of the surface can roughly be defined as the time it takes the surface to move a wavelength — 
in the present case X/v. This thought experiment is useful to understand the limitation 
in coherent integration time for a moving surface like the ocean. An advantage of using 



23 



a synthetic wavelength using two GNSS frequencies is that this coherent integration time 
should be longer. We will verify in a moment that this is the case. 

An associated concept to the coherence function is the structure function of the phase. 
This is defined by 

s (T) = m + T)- ( p(t)) 2 )f (6-3) 

It gives us a measure of how tied up are two points in the phase as we vary their temporal 
separation. If this phase drift can be approximated by a random walk stochastic process 
[ Herring, 1990^ , 

s{T)=d-T (6.4) 

we can use the associated "Geophysically Induced Phase drift" (GIP) rate \fd as a geo- 
physical parameter if it correlates well to U10 in our models. This parameters has units of 
cycle per square root second. We will check this below (but let us anticipate that it does). 
The coherence time and the structure function are related. In fact, 

(C(T) 2 ) = -L j T f 1 \jW>-Wy))dbd£ '. (6.5) 
T z Jo Jo 



If we assume Gaussian statistics for the phase fluctuations \ Thomson et al, 199L ] we can 
simplify things considerably, 

(C(T) 2 ) = -L [ T ( T e-'V-Wdtdtf. (6.6) 

1 Jo Jo 

2 rr 'l-l=)e-<rWdr. (6.7) 



Tin \ T 



For a random walk process, then, 

(C(Tf) = ^J o T (l-f)e- d ^dT. (6.8) 



See Figure I6J for a graph of this function for various drift rates. 

Note that in this model (C(T) 2 ) decreases to zero as time increases. This fact alerts us 
that there is something amiss: we see in our simulations large phase excursions (certainly 
much larger than a cycle) , but there is coherence left at large times — see Figures |6.2|, pi 



and 6.5. How is this possible? The Gaussian model for phase fluctuations does not really 
apply to our situation, however. Gaussian means that phase excursions probabilities are 
bell-shaped, of course. Are they? In our simulations the phase of the reflected signal gathers 
about a point in the complex plane — the field does seem to have a Gaussian distribution 
about an average point. The cumulative phase does not have to have such a distribution, 
however. In fact, it could conceivably indulge in arbitrarily large excursions while spending 
more time in that average field point: this would still result in a Gaussian field distribution. 
However, the probability distribution for such a phase history would not be a Gaussian 
distribution (centered about 0, say). It would consist of a series of Gaussian-like humps 
spaced one cycle apart, with the one about zero the largest. In this situation we could 
well have a very long coherence time (since the phase modulo 2tt could be a rather narrow 
Gaussian distribution about a point) while the structure function may show high drift. This 
type of behavior is evident in our simulations when plotting the field points (see Figure |6l 



Note that the interferometric signal has a stronger average field value. See Figure 6.7 for a 
histogram of the L25 cumulative phase illustrating the multiple-hump distribution that we 
alluded to. 
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<C>2(tau)> with a Random walk, s(t) = d tau. sqrt(d) = 0.5, 1.0, 2.0, 3.0 




0.0 I , , , I , , , I , , , I , , , L 

2 4 6 8 

Time, in seconds. 



Figure 6.1: (C(T) 2 } for a random walk with drift rates of \fd =0.5, 1.0, 2.0 and 3.0 cycles 
per square root second. 
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Critt mta= i- JBIUTJ p*r 





Crift rat?- £.1 pggg pir rrthw?; 




Tht\ [ton*] 
Drift ratal = j Jfl8ft7 ^Hm par Ecrtfratj 





Figure 6.2: This is a plot with the structure function (left) and the coherence function 
for the case U10=6 m/s, 43 m side 10 cm resolution simulation at 8 m height. This is 0.4 
meters height standard deviation. From top to bottom, LI, L2, L5 and L12, L25 fields. 
The drift rates are 4.2, 2.1, 1.1, 6.2, 2.5 cycles per square root second. Coherence times 
0.04, 0.06, 0.08, 0.09 and 0.3 seconds. 
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Figure 6.3: This is a plot with the structure function (left) and the coherence function for 
the case U10=l m/s. From top to bottom, LI, L2, L5 and L12, L25 fields. The drift rates 
are 1.9, 0.9, 0.8, 1.4, 1.3 cycles per square root second. Simulation size is 19 m side, 10 cm 
resolution, and height standard deviation is 3 cm. 
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Figure 6.4: This is a plot with the structure function (left) and the coherence function for 
the case U10=3 m/s. From top to bottom, LI, L2, L5 and L12, L25 fields. The drift rates 
are 0.8, 1.0, 1.0, 1.0, 1.5 cycles per square root second. Simulation size is 19 m side, 10 cm 
resolution, and height standard deviation is 8 cm, receiver height 8 m. 
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Figure 6.5: This is a plot with the structure function (left) and the coherence function for 
the case U10=4 m/s. From top to bottom, LI, L2, L5 and L12, L25 fields. The drift rates 
are 3.7, 2.1, 1.6, 2.1, 1.8 cycles per square root second. Simulation size is 19 m side, 10 cm 
resolution, and height standard deviation is 12 cm, receiver height 8 m. See next Figure 
for further analysis of this particular simulation. 
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Figure 6.6: Refer to the previous figure. This is a plot of the field in the complex plane 
for Li2, L5 and L25, and it serves to illuminate the difference between the structure function 
and the coherence function. The field mean for L2 is (-0.06,-0.1), for L5 is (-0.1,-0.1), and 
for L25 (-0.6,-0.20), much larger. The standard deviation is about 1 for all of them. In this 
19 meter simulation, U10=4 m/s, height standard deviation was 12 cm, receiver height 8 
m. 
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Figure 6.7: In these two plots we depict the phase and phase histogram for the L25 field 
in the previous figure. The highest peak corresponds to the initial zero phase. 
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Notice that the above discussion is relevant to understand the coherence properties that 
were found in the Bridge Experiment reflected signal (integration times where of 10 ms, 
short enough for the static approximation to work in that geometry). In longer time scales 
the accumulated phase will pick up a random component, but if there is an overall phase 
change (due for instance, to a tide) it should be distinguishable from the random part 
(which should not add nor subtract to the overall phase on the average) . 

6.2 Doppler spread of the reflected signal 

Another thought experiment is useful. Imagine that the surface height is varying as a 
function of time as 

((t) = at + Ar(t), 

where a is a constant and the second term adds a random, noise-like displacement between 
—A and A. We can think of this as the effect of a tide superimposed on the wave mo- 
tion on the sea. The scattered phase will not be coherent if A is large enough, but the 
accumulated phase will certainly contain information about a. What happens if there are 
many scatterers, how fast is the residual phase varying? This will give us a feel for the 
Phase Lock Loop (PLL) bandwidth required for tracking. It is not hard to see that the 
bandwidth required will be less than 100 Hz (assuming that the sea surface vertical speed 
is not greater than 10 m/s), in a static geometry (or if a static- WAF is used). The largest 
frequency present in the signal cannot be greater than the contribution due to the fastest 
moving surface patch — the result is always the linear combination of such contributions, 
U(t) ~ re*^w = Y^=i e% ^ j ^ ■ Now, if the maximal orbital velocity for a wave is 10 m/s 
(see Figure [T^, this implies a Doppler of 2 x 10/ A ~ 100 Hz (remember that the signal is 
bouncing off the surface). Thus, a sampling rate of at least 200 Hz is required to retrieve 
the phase of the reflected signal (even if the SNR is high enough). In our simulations we 
have seen that it is safe to sample the signal at such rate. 

Here is another useful question: given that a receiving LEO is moving at, say, 7 km/s, 
would it be possible to integrate coherently (correlate) the direct and reflected signals for 
1 ms, if the LEO antenna has a footprint of 10 meters? Note that this "eye on the sea" is 
moving very fast over the surface, and this is introducing a random, time-dependent phase 
to the signal. The answer is no. But if the footprint is enlarged enough, yes. This may 
be an important problem for the PIP, but there are ways through which it could be fixed. 
One obvious one is to use a WAF that sticks to a given surface patch, or that it moves 
slowly over the sea surface. Some of these ideas are discussed within the context of the 
GNSS-OPPSCAT project (Report WP3320). To get an idea of the induced Doppler, we 
can reason that the scattered field is just a sum of the type 

N 

F = R(t)e 1 ^ «5>^)e*W, 
i=i 

which already looks like a Fourier sum if we approximate ri (t) as constant and the exponent 
by its first derivative. Thus, we can guess that the maximum frequency in the spectrum is 
simply Doppler due to the fastest moving patch on the ocean. As an upper bound we can 
use 10 m/s — see Figure [6~9| . Multiplying times 2 and dividing by A we obtain a scattering 
Doppler of about 100 Hz. 
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To be a bit more precise, let us write 

F = f e iq * z d 2 x = [ e iq ^ Zi,+zt) d 2 x, 
■> Awaf JA W af 

using the simpler backscattering expression in the Fraunhofer zone. Notice the restriction 
of the integral to the (static) WAF zone. Let us now rewrite the transform of the field: 

F{uj) = J Fe iujt dt = J e iq * z ° J e^'^dtdx 
~ / e- iqzZ °d 2 x, 

J A(qz=u>) 

if we assume that vertical motions take over a "long enough" time compared to the time 
scale associated to the frequency we are interested in. This is saying that the spectral 
component at a given frequency is proportional to the area on the surface moving at the 
right speed to produce that frequency. 

We can extend this reasoning to the moving-receiver case. There are two limiting 
situations we can think about. One is the "high gain" situation. Let us first discuss this 
situation — we will get back to the low gain situation at the end. Assume for now the 
ocean is frozen. In this situation the WAF zone is very small (due to antenna or processing 
gain) and coherence time is definitely influenced by the fact that as the receiver moves the 
contributing surface is changing. If the WAF zone is smaller than the correlation length of 
the ocean (certainly an extreme case, but never mind), the coherence time is limited by l/v, 
where v is the receiver speed and I is the correlation length of the surface — assuming large 
enough height deviations. To be more precise, if we use a replica with fixed Doppler and 
Delay, the "eye on the sea" will move over the ocean at a speed similar to the receiver's. 
This will introduce additional bandwidth in the carrier if the ocean is rough — even if is 
frozen. This bandwidth will be mainly proportional to ocean surface rugosity, not to ocean 
surface motion. It can be very large (kHz), depending on the receiver speed and the WAF 
area size. Following the above reasoning, and keeping in mind that the WAF zone is now 
a moving filter on the ocean surface, 



F(t) = f e iq * z Wd 2 x= [~ WAF (x-vt)e iq * z Wd 2 

_ J cWAF^j e iq z z(x+vt) J. 



#X 



= f e iq ^ +i!t) d 2 x 

Now let us look at F(t): 

F(t) = v- J V U<i^+vt)\ d 2 x = iJ . I e iq z z(*+vt) ft dl ( ) 
JA WAF V ' JdA WAF 

where the last is an integral over the boundary of the WAF area and n is the normal to 
the boundary. The result illustrates the "edge effect", i.e., that the change is just due 
to the change at the boundary of the WAF area. We have used a 2D version of Green's 
theorem in the last step. This equation says that the change in the field is proportional to 
the variation of the field contribution at the edges of the WAF area (properly mapped in 
the velocity direction). The faster the velocity and the larger the difference, the larger the 
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rate of change of the field. To relate this to the spectrum, note that this will contribute 
a high frequency componene to the spectrum. Whether this is a relatively large or small 
contribution depends on the total field, which is proportional to the total area. Thus, 
for large WAF zones this high frequency effect contributes a small portion of the total 
spectrum. 

If we allow for a moving surface, z = z(x, t) it is readily seen that the complete result 
is the sum of two distinct effects: 

F{t) = v- [ e iizz(x+vt,t) fidl+ f e waf (x - vt) d t d 2 x. (6.10) 

JdA WAF J v ' 

Let us now return to 

F(t)= [ e iq ^d 2 x= ( E WAF {x-vt)e iq * z ^d 2 x (6.11) 

JA WAF (t) J 

and compute the Fourier transform: 

F{u) = I r dt e iujt E WAF (x-vt) e iq ^d 2 x. (6.12) 

J J — CO 

To evaluate this, let us assume for simplicity and without loss of generality that E WAF is 
a box of size S moving along in the x direction with velocity v x : 



Then 



Hence, 



E WAF (x-vt) = E x (x-v x t)-E y (y). (6.13) 

/oo c („.\ p iu>x/v x 
dt e iujt E WAF (x - vt) = 2^^- sm(u;S/v x ). (6.14) 
-oo W 

F(u) = 2sin (° jS / v ^ f Ey{y) e Wtfc e i ^d 2 x. (6.15) 

The characteristic time is roughly given by S/v x , as can be seen from the multiplying sine. 
Frequencies higher than v x /S are supressed by this factor. The power in frequencies smaller 
than this are modulated by the horizontal roughness of expiq z z(x) in the scale defined by 
v x /lo in the direction of motion: if the surface does not vary in that scale, there will be 
little power at the frequency v x /S. 

In this approximation we have implicitly assumed a very small WAF in relation to the 
geometry, since we have worked in the Fraunhofer approximation. This means that all the 
points in the reflecting patch have the same geometry and generate the same geometrically 
induced Doppler. In recent experimental conditions this is hardly the case, and the Doppler 
of the reflected signal is dominated by the Doppler amplitude of the WAF and glistening 
zones. If the glistening zone is large enough, then radiation will be received from the entire 
WAF. The Doppler span of the WAF zone determines then the Doppler spread of the 
received field. This is the basic mechanism for Delay-Doppler mapping and the basis for 
the relationship betweeen waveforms and the characteristics of the sea surface. Even if we 
restrict the signal to the first chip in the C/A code, there will be quite a bit of geometric 
Doppler. This second source of Doppler is in principle removable by many means: high 
antenna gain or processing gain. By any of these means we can reduce the patch to the 
Fraunhofer zone mantaining SNR. It will then become increasingly important to use SAR 
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T, (ras) 
-*- % \ ' ' vkj ) 


Ocean motion (Hz) 


Geometric (Hz) 


Edge (Hz) 


i 


200 


1000 


1 


5 


200 
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Table 6.1: Back of the envelope calculations for Doppler spread induced by ocean motion, 
geometry and edge effects. Note that above 20 ms edge effects force the use of focusing. 



techniques to focus and eliminate the edge effect we just discussed. We plan to revisit these 
calculations extending them to the Fresnel zone in future work. 

The ultimate limit to the Doppler width and therefore to the coherent integration time 
of the signal is ocean motion-that cannot be anticipated and compensated for. 

The size of the first Fresnel zone for a receiver at 350 km is about 2^/2h\ or about 1 km. 
In this narrow area, the geometric Doppler spans about 200 Hz edge to edge (the receiver 
looking down and is moving along at 7 km/s) and the corresponding integration time is 5 
ms. This Doppler is of the same order of magnitude of ocean induced spread. Integration 
times longer than 10 ms will certainly begin to be sensitive to ocean motion. The edge 
effect we discussed earlier will also start to play a role at longer integration times : for a 1 
km WAF (roughly corresponding to 5 ms integration time) the corresponding frequency is 
of 7 Hz — still too small to make spotlight processing necessary. This may be the best way 
to retrieve ocean induced Doppler spread. See Table- 3.1 for a summary of these ideas. 

We have two different new ways to extract sea-state information from the Carrier band- 
width, by playing with the eye location. One the one hand, we can use a static WAF to 
measure vertical velocities — this will entail using SAR techniques from moving platforms. 
On the other, we can measure surface roughness by letting the WAF drift over the surface 
at a speed of choice. This is just a matter of tuning the matched filter appropriately, and 
it could even be done from a static platform. 
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Figure 6.8: The ocean surface and vertical velocity for a wind of UW = 8 m/s, Elfouhaily 
spectrum. The contour denotes the zero value in each case. 
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Figure 6.9: Ocean heights and vertical velocities using the Elfouhaily et al. spectrum for 
wind speeds of U10=l,3,5,..,17 m/s. Note that due to limitations in the simulation size 
values for U10 above 12 m/s suffer from saturation. 
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Chapter 7 



Statisitical Properties of the 
Reflected Fields 



7.1 Analysis of field correlations 

In this section we address the question of the correlation between the fields at different 
frequencies (say L\ and Li and the future L5) for different sea conditions. An approach to 
this problem is to rewrite the Fresnel integral using the "zone" concept. That is, we classify 
areas in the surface according to their distance to the receiver. We assume we are in the 
Fraunhofer-emitter zone, so the distance to the emitter plays no role in this discussion. The 
idea is then to rewrite the field as a sum of field contributions from equal delay (r)zones: 

e iqr A'(r)dr. 

■ min 

This idea can be found in [Berry, 197i\. The function A'(r) is an "areal" density function 



that takes into account how much area contributes to each zone. It's exact expression is 
not of immediate concern. In the case of a flat surface it just becomes one. In order to 
extend the integration over the whole r-axis, let us simply extend the definition of A'(r) to 
be zero outside r m i n and r max . We then have, 



/oo 
e iqr A'{r)dr. 
-00 



We can now read: the field is the Fourier transform of A'{r). And the question about the 
relationship of the field between frequencies becomes a question on the relationship between 
the Fourier components of this areal function at different frequencies. 

Let us forget for now the meaning of the areal function, and just think in terms of the 
Fourier transform of a function with limited support, since we are looking at the Fourier 
transform of a function with support L = r m i n — r max . The larger the support of the 
function, the more small-scale structure can be found in the Fourier transform — everything 
else the same. We are accostumed to thinking about this fact in the opposite domain: 
a function with a lot of small time-scale structure will have a large bandwidth — a large 
support in the frequncy domain. 

Intuitively, the larger the support two given nearby frequency components will end up 
sampling the pulse at more separated regions. It is useful to think about two sine functions 
of slightly different frequencies, running side by side. As a rule of thumb, we expect that the 
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correlation between these Fourier components will be sensitive to the correlation function 
of the (areal) function at a distance (A A) * L/X mean . This is because A A is the sampling 
distance difference gained per cycle, and *L/X mean is the number of cycles available in a 
region of size L. Therefore, before carrying out any calculation we can say that if the 
areal function is correlated at the distance dictated by this (maximal possible) frequency 
separation, then we will see correlation in the fields. It helps to look at the separation 
between the equi-delay contours in the surface of integration. If this separation is smaller 
than the correlation length of the surface, then high correlation between the fields is possible. 
There is an extra consideration, however. We have been assuming that the separation 
between contours is constant, but this is not the case. That is, we have been implicitly 
assuming that the contours are fixed on the surface, and all we have considered is the 
correlation between this contours. But in fact, as the ocean moves, the contours move, 
especially with large seas. This introduces an additional source of decorrelation, which gets 
worse with larger seas. The time-correlation between the fields is given by 



(A' 1 (r)A' 2 (r')) t e iqir - iq2r ' drdr 1 '. 



This is a difficult beast to deal with, so we will change strategy. 

The correlation between the fields can also be calculated in another manner. It will be 
useful to approximate things entirely in the Fraunhofer zone, i.e., with both transmitter 
and receiver in the far field, assuming Gaussian statistics. The field in the nadir case is 
given, in 1-D and up to a constant, by 



U = q e iqz{x) dx. 
J A 



Now, 



and 



|c/p 



AJA 



e iQ(z{x)-z(x')) dxdx l^ 

( e «9l*W-*l* »)dxdx' ', 



I AJA 

2 I I /Jq(z(x)-z(x')) 



(7.1) 

(7.2) 
(7.3) 



Now, if we assume Gaussian statistics, with 

1 



P(z(x) — z{x )) 



27R7V 1 - P 2 ( X 



: CXp 



z(x) 2 — 2p(x — x')z{x)z(x') + z 2 (x') 
2<r 2 (l- p 2 (x-x')) 



(7.4) 

we can carry out these calculations explicitly. Here a is the height standard deviation 
from the (zero) mean, and p is the correlation function of the surface. It is healthy to 
keep in mind that in the gaussian case, there are the only parameters, together with the 
frequencies, that can appear in the final expressions. If we further assume a Gaussian 
correlation function with correlation length I, the answer to any question must be expressed 
in terms of dimensionally meaningful expressions contantaining, a, I, Ai and A2. 
For the case at hand, the key result is ( ^Beckmann et at, 1963^ , p. 190) 



fJq i z{x)-\-q2z{x')\ 

This implies, for instance, 



exp 



-a (q 1 + 2p(x - x)qiq 2 + q 2 ] 



Q1Q2 / exp 
I A J A 



-a {q x - 2p(x - x')qiq 2 + q 2 ) 



dxdx' . 



(7.5) 



(7.6) 
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In the following we use a simple approximation, 

p{u) = 1 — \u\/l, \u\ < I, 



(7.7) 



else zero. In this simple case we model a Gaussian correlation function by a linear approx- 
imation. For simplicity we will also work in the 1-D case. The first step now is to define 
u = x — x', and v = x + x'. The Jacobian of this change of variables is 1/2. The integral 
becomes (let A be a 1-D area from — L to L) 



V2L rV2L-u ^udv 



Q1Q2 



/V £1-' 
-V2L J- 



V2L Ql 



-\/2L J-V2L+U 2 
r\/2L 



exp 



-^°" 2 (9i - 2p(u)q 1 q 2 + ql) 



Q'2 



V2L 



du exp 



-^a 2 (ql + 2p{u)qiq 2 + ql) 



(7.8) 
(7.9) 



Now, using our simplified correlation function model (and we will henceforth assume that 
L > I), we find 

(UiUZ) = 2V2L (Jj (l _ e - 2 n^ e -aW/2 + qiq2 (^2 L _ i^ 2 (g( +q ^)/2^ _ (7 _ 1Q) 

For L-band, q = 2k is about 60 per meter. If a is greater than 0.1 m, a 2 * q 2 is greater 
than 40, and all the exponential terms in this expression vanish, except for exp[— a 2 Aq 2 /2]. 
That is, for a > 0.1 m, 

2V2LI _ a 2 Aq 2 /2 



It also follows that 



(UiU* 2 ) 
C(U 1 ,U 2 ) = 



^/(UiU?){U 2 UZ) 



~ e 



-a 2 Aq 2 /2 



(7.11) 
(7.12) 



1 .2 
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Figure 7.1: The coherency factor C(Ui, U 2 ) (dashed) and C(U 2 , U5) for the case / = 5 and 
L = 20. There is not much sensitivity in either case to L or /. 
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According to this result, the critical parameter in any reasonable ocean state is simply 
the significant wave height. Imagine for instance an ocean with a very large correlation 
length and large height standard deviation. This is mirror-like ocean with a global up- 
down displacement. It is quite clear that ih this example (t/it/J) will be zero, since this 
product will be a number in the complex plane with a reasonable magnitude and a rather 
random phase. The average of such a complex number is zero. 

The general result, relevant for smaller significant wave heights, is 



C(U 1 ,U 2 



/ 

171 



1 



+ q\{y/2L - l)e 



1 



Hi 



+ q%(V2L-l)e 

(7.13) 

We show a plot for the case I = 5 and L = 20 in Figure 

Although we did not discuss it, it is rather immediate that the average field for the pure 
fields are governed again by exp [—a 2 q 2 /2] — this can be read off Equation 7.5 assuming a 
zero correlation function to infer (exp [iqz(x)\) = exp [—a 2 q 2 /2]. 

We have checked this trend in our data. In general, the interferometric combination 
average fields are larger. For instance, for a cr^= 79 cm simulation (50 seconds at 0.0025 
temporal resolution, 44 m size, spatial resolution 10cm, U10=1100cmpers) — see Figure 7.4, 
the LI, L2, L5, L12, L25 average fields were: 



LI 
L2 
L5 
L12 
L25 



(+0.03,-0.01) 
(+0.00, -0.06) 
(-0.04, +0.03) 
(-0.02, +0.09) 
(-0.12, +0.07) 



while for U10=4 ms, a^= 12 cm simulation, 19 m side (see Figures 
more details about this simulation) 



.5 , 3.6 and 



for 



LI 
L2 
L5 
L12 
L25 



(+0.05,-0.03) 
(-0.06,-0.12) 
(-0.04, +0.10) 
(-0.09, +0.12) 
(-0.63, +0.20) 



This is where the real strength of PIP becomes evident. Moreover, low-pass filtering (see 
next section) does not change these numbers much — the average field is robust. 

The Fraunhofer approximation has given us a useful means to understand the correlation 
properties of the signals in different frequencies. We should note, however, that in the GPS 
case we are not in the Fraunhofer zone — we are in the WAF zone, which is is usually 
substantially larger. Understanding all the implications of the extension to the WAF zone 



involves the analysis of the "beast" in Equation which is substantially more complicated 
by the factors inthe areal function. This is left for future work. 



7.2 Filtering 

Another way to understand and take advantadge of the coherence proporties is to focus 
on the low frequency aspects of the reflected signals for the static case. As we mentioned 
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above, in the static situation the coherence function of the scattered phase is related to the 
zeroth Fourier component of the field — see Equation 6.1 . Filtering the signal to retrieve slow 
varying geophysical signals is thus a plausible approach. As an example, we show in Figure 
7.3 the fields with and without filtering at 0.5 Hz (relevant for the case of a LEO). The 
result is quite good, as can be seen. This behavior is also expected from the analysis above 
of the average field above, of course. The interferometric combination is again superior, 
especially L25. For the satellite of aircraft case, good modeling of the receiver position and 
filtering can be used to remove all but the slowly varying geophysical signals from the phase 
drift. 

Important points: filtering doesn not seem to affect the field mean value (a very desirable 
result). It does significantly reduce the scatter, however. For instance, in Figure |7^, we 
have the fields for a 3 m/s wind. The standard deviations of the field go from 0.87, 0.82, 
0.79, 0.73, 0.68 to 0.51, 0.51, 0.49, 0.48, 0.29 (for LI, L2, L5, L12, L25) after 0.5 Hz filtering 



of the fields. See also the results for higher wind speed is Figures [7.3| and 7.4. The mean 
field remains virtually unchanged after this process. 
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Figure 7.2: These are plots of the fields and phases for the case U10=3m/s (cr^ = 18 cm), 
with an ocean patch of 43 m, resolution 10 cm. The top is the original time series, the 
second is with a field filter at 0.5 Hz. Temporal resolution in all the simulations is of 0.0025 
seconds. 
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Figure 7.3: These are plots of the fields and phases for the case U10=6 m/s (cr^ = 42 cm), 
with an ocean patch of 43 m, resolution 10 cm. The top is the original time series, the 
second is with a field filter at 0.5 Hz. 



44 



LI Dh 




L2 pHq 



D ID 2D 33 +0 3D 



Magnitude L1 




D 10 2D S3 40 ED 
Time, BKDnta 



FFTCfldd) 




Hi 



ptlQEH 




Magnitude L1 




D 10 2D S3 40 ED 
Time, BKDnta 

FFTCfldd) 




Hi 




D 1'3 H' M HI J[ 



Magnitude- L2 




D 10 20 S3 40 ED 



FFTtfMd) 




1O3J0J34O50 

Hi 




D 10 20 30 -ID 3D 
Magnitude L2 




LI 2 ph38 



D 1G ZD 3D W 30 



Magnitude L5 




D 13 20 3D 40 ED 
Uma. ucarHta 



FFT{fl9ld> 




1C2M04O5O 

Hi 




D 1Q 20 3D 40 SO 
Tims, nconda 

Magnitude L5 



D 10 20 S3 40 ED 
Hmu. iKDnA 



FFTtflgld} 



D 13 20 SD 40 ED 
Uma. ucarHta 



FFTtflBld> 






13 2D 3D -ID 30 



MagnitLido L12 




D 13 20 SD 40 EC 



FFTfflsId) 




1O3t'304O5t' 

Hi 



LI 2 phow 




1 0.6 




|,.D 




i 2J3 
U 

a 






^ a& 




js 1JD 


f D.D 


i ai 

f ttB 




i oj 


lira: 




f OJO 



13 2D 3D 40 30 
Magnitude L12 



D 13 20 3D 40 EC 
ITniu. 3tKandi 



FFT(fldd) 




L25 pha 



10 





















f 


-10. 






-1B 






-2D 







13 20 3D 40 SO 
Urns, ¥KCfidi 



Magnituda L25 




1jO 
fa O-fl 

T 

a 

02 
M 



D 13 20 SD 40 ED 
limn, ucwith 

FFT(fleld) 



10»3S4O5O 

Hi 




13 2D 3D +0 BO 
TTma. ¥Q«fidi 

Magnitude L25 





2J0 








£ 


1-5 










1JD 


u 




OlB 




i 


OjO 





1J0 
Ofl 
fa 0_fl 

T 

a 

02 
M 



D 13 20 ID 40 ED 
limn, ucwith 

FFT(fleld) 



1O3JOJ34C50 

Hi 



1O2M04O5O 

Hi 



10310304050 

Hi 



1OW304O5O 

Hi 



Figure 7.4: These are plots of the fields and phases for the case U10=llm/s (<j£ = 79 
cm), with an ocean patch of 43 m, resolution 10 cm. The top is the original time series, the 
second is with a field filter at 0.5 Hz. See the next figure for the corresponding (unfiltered) 
coherence and strucuture functions. This is complete simulation, with the ocean going from 
flat to moving to flat. Observe the picked up winding number. Filtering at 0.5 Hz did not 
entirely cure the problem in this case. 
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Chapter 8 



Simulation experiments with 
winding number and GIP 



Consider the following scenario. The receiver is static. Initially the ocean is calm, totally 
flat. Then a disturbance gradually appears, peaks, and slowly disappears. The ocean final 
state is the same as the original state: total flatness. We have run several simulations with 
this scenario. The result is that the winding number of the field history is non-zero. That 
is, the field can loop several times arond zero. 

The first example is in Figure 8.1 and accompanying Figure 3.2. How is this possible? 



As can be seen in the example in Figure |8.3| , the sum of zero-winding-number curves can 
result in a winding number 1 curve. Since the resulting total field is a linear combination 
of such fields, all that we need to show is how the moving ocean can generate such winding 
number zero curves. This is easy: the only requirements are that a) the initial and final 
phase be the same, with the phase moving in between, and b) that the field magnitude 
increase or decrease before returning to its original value. It is possible to imagine examples 
based on GO reasoning (a usefull approximation, at any rate) to obtain this behaviour. 
Suppose, for instance, that the distance from a specular point to the receiver decreases and 
then increases, while the radius of curvature decreases then increases. It is easy to see that 
this yields a clockwise loop in phasor space (with the usual angular convention used in the 
complex plane). 

We have seen the phase dispersion effect, very much related to winding number (which 
is a special case of the first when initial and final states are the same), in all wind conditions, 
with U10=2 m/s and up. At first sight this poses a significant problem in the use of phase 
for altimetry. How can we use a "denting" measure for altimetry? It is now very important 
to characterize the statistical and geophysical properties of this phase drift (Geophysically 
Induced Phase drift, or GIP for short). 

8.1 The uses of winding number 

Winding number, or more generally, phase dispersion under a moving ocean is certainly 
related to sea state characteristics. An interesting GedankenQ imagine a satellite tracking 
a reflection over a smooth region, then rough, then smooth again. The winding number 
induced will be related to the roughness encountered. We have already seen that, in general, 

1 Gedanken (thought) experiments are little mental exercises for exploring hypothetical experiments. The 
terminology comes from Einstein's work, who is said to have engaged in them even when he was very young. 
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Figure 8.1: Simulation: 20 meters side, res=10cm, max ocean 4 m/s U10 with a peak 
of 12 cm, sine fourth modulation. Note that in the bottom row we present the FFT power 
of the normalized field. The zero component is basically the coherence function mentioned 
earlier. The interferometric component shows a higher degree of coherence. 



the phase drifts without a specific chirality in our simulations. It remains to be seen if non- 
gaussian effects in the surface can lead to chirality, or preferred rotation direction. If this 
is the case, this phenomenon could conceivably be used for geophysical measurements. We 
have already seen that GIP is related to wind speed in our simulations. 



47 
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Figure 8.2: In this Figure we illustrate the L2 field phasor in Figure 8.1, which has a 
winding number of -1. 




Figure 8.3: Winding number magic: four curves with winding number yield a winding 
number 1 curve after being added. The trick can be done with 2 curves. In this case they 
all turn in the same direction. Changing that will result in an ellipse. 
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Chapter 9 



Recommendations for 
Post-processing Procedures 

Based on the work carried out so far, we are prepared to outline some of the possible 
incarnations of the PIP concept. What we must emphasize is that the strength lies in the 
noise cancellation in the interferometric combination. Filtering is best carried out after this 
combination. 

9.1 The PIP instrument: architecture and data products 

Based on our results, we have the following suggestions. 

1. Use two wideband (200 Hz for the static receiver, possibly more for the moving case) 
PLL's to extract the phase in each frequency. 

2. Combine the phases and low-pass filter to recover altimetric trends. The normalized 
interferometric field is more coherent, so it can be used to track slow- varying altimetric 
changes (such as a tide) . We have seen that the interferometric fiels is a more robust 
source, with a better defined phase. 

3. Use phase drift as a geophysical parameter More generally, the phase and power time 
series are probably rich information source of sea state. 

To be more precise, the following recipes for PIP implementation can be conceived: One, 
the most natural one perhaps, is the Single PLL PIP (SPIP) 

1. Multiply the two signals, form the interferometric signal. 

2. Filter from to 0.5 Hz 

3. Use a single PLL to extract the phase. 

The disadvantadge is that the field magnitude will oscillate quite a bit. This may not be 
important, given that a filter will be used. An alternative would be to used two PLLs 
(MPIP) 

1. Extract the phase in LI and in L2 

2. Combine the two phases 
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3. Filter from to 0.5 Hz. 



This is probably similar to the previous case, but it takes the field magnitude out of the 
picture from the very beginning — probably not a good thing. Finally, a non-contender is 
to 

1. Filter from to 0.5 Hz 

2. Extract the phase in LI and in L2 

3. Combine the two phases 

This is a bad option, it makes no use of the jitter cancellation across frequencies. 

We have focused on the static case. For other situations we just need to work with 
excess phase and carry out the same filtering process. 
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Chapter 10 

Conclusions and Future Work 



10.1 What we have found: PIP's superior performance 

1. The mean field is largest for the PIP case L25 for all sea conditions. In general the 
mean field is proportional to exp [— a 2 q 2 /2], where q = 4ir/\ and A is the pure or 
interferometric wavelength. 

2. The coherence integration time is consequently substantially larger for the interfer- 
ometric combination, especially L25 (if the mean sea level scatter is very large the 
effect disappears, of course). The interferometric phase is coherent in a reasonable 
range of sea conditions, unlike the single frequency phases. This implies that the PIP 
interferometric combination is superior to extract altimetric low frequency trends in 
the phase, as we discussed at the beginning. 

3. The phase in our simulations behaves like a random walk, and the drift rate is directly 
related to sea state: in general, the larger the wind speed, the higher the drift rate. 
This implies that a system capable of tracking the phase of the reflected signals can 
provide, aside from altimetric measurements, sea state information. 

4. Filtering the received field does not affect its mean value, while decreasing its scatter. 
It can probably be used for altimetric purposes, as it removes noise but leaves the 
slow- varying geophysical signals alone. 

5. No long-term trends in phase drift chirality have been detected in our simulations 
with Gaussian ocean models based on ^Elfouhaily et at, 199% . We have performed 
simulations up to 2 minutes long. This is good news for our altimetric efforts. It seems 
that despite the phase drift present, and the winding number it can lead to, in practice 
this effect should not introduce systematic effects on altimetric determination. It is 
important to check that the simulations are correctly imitated by nature, of course. 
Among other things, non-Gaussian effects in the sea could conceivably lead to real 
drifting. 

The single most important parameter in the correlation between the fields at different 
frequencies is the significant wave height. This has to be compared to the synthetic wave- 
length. For rough seas, the correlation bethween the fields (and therefore interferometric 
coherence) in different frequencies disappears, and the coherence time goes to zero even for 
the interferometric combinations. In calmer ocean conditions, however, our results indicate 
that the interferometric combination remains coherent. 
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10.2 What we'd now like to know: future work 



Due to the limited scope of this study we were not able to cover all the interesting aspects 
of this problem. Future theoretical work includes the extension of the present results to the 
case in which the receiver is far away from the ground and is moving. We would also like 
to understand the impact of changing the ocean statistics from gaussian to non-gaussian. 
Another important aspect that deserves further research is to extend the analysis to a 
bistatic situation. So far, our simulations have only considered the monostatic case. We 
intend to continue this work within the scope of PARIS-a: 

1. Examine more realistic situations (higher, faster). This includes, in particular, un- 
derstanding the issues that will arise in the aircraft and LEO scenario, and it will 
require larger simulations — in fact it may impossible to simulate the ocean with the 
required size and resolution, as the WAF zone are increases linearly with height. 

2. Understand the theoretical issues involving the correlations and Doppler spread of 
the reflected fields without using the Fruanhofer approximation for the receiver. 

3. Extend the analysis to the bistatic situation. 

4. Carry out a detailed study of phase drift versus wind speed. 

5. Look deeper at the advantages (which we have seen and are expected) of using even 
closer frequencies. This may result in recomendations for GALILEO, for instance. 

6. Understand the impact of non-Gaussian ocean effects. 

7. Understand better the possible sources of chirality in GID. This is important to assure 
altimetric accuracy. 

8. Obtain and analyze data for the static situation. 

9. Understand the effectiveness of GO or its higher order corrections. 

It is rather clear that experimental work is needed in this field. The most important 
experiment for the concept at this point, is a redo of the bridge experiment: a static 
receiver over the moving ocean, using both Li and L2. The desired product from such 
an experiment from the point of view of the present work, would time-series' of the fields. 
MEATEX campaigns will also provide aircraft data which will be very useful to understand 
some of these issues. 



52 



Chapter 11 

Acknowlegments 



This is a document produced for ESA under ESTEC Contract No. 14071/99/NL/MM. The 
author's are grateful to ESA for funding and permission to freely distribute this document, 
and are especially grateful to Manuel Martin-Neira of ESTEC, the technical officer in charge 
of this contract, for very valuable comments and suggestions during the course of this 
research. . 



53 



Bibliography 



[Auber et al., 1994] Auber, J.-C, Bilbaut, A., Rigal,J.-M., Characterization of multipath 
on land and sea at GPS frequencies, Proceedings of the 7th International Technical 
Meeting of the Satellite Division of the Institute of Navigation, part 2, ION-GPS-94, 
Sept 94, pp. 1155-1171 

[Beckmann et al., 1963] Beckmann, P., Spizzichino, A., The scattering of Electromagnetic 
Waves from Rough Surfaces, Artech House, Inc., Nordwood, MA, 1987. 

[Berry, 1972] Berry, M.V., On deducing the form of surfaces from their diffracted echoes, 
J. Phys. A: Gen Phys 5, 272-91. 

[Born & Wolf 1993] Born, M., Wolf, E., Principles of Optics, Sixth Edition, Pergamon 
Press, 1993. 

[Caparrini, 1998] Caparrini, M., Using reflected GNSS signals to estimate surface features 
over wide ocean areas. 

[Chapron et al., 1999] Chapron, B., Kerbaol, V., Vandemark, D., Elfouhaily, T., Impor- 
tance of peakedness in sea surface slope measurements and application, submitted to 
Journal of Geophysical Res. , March 1999. 

[Cox and Munk, 1954] Cox and W. Munk, Measurement of the roughness of the sea surface 
from photographs of the Sun's glitter, Journal of the Optical Society of America, vol. 
44, pp. 838-850, Nov. 1954. 

[Daout et al, 1999] Daout, F., Schmitt, Characterization of the Bistatic scattered distri- 
bution, 1999. 

[Elfouhaily et al, 1991] Elfouhaily, T., Chapron, B., Katsaros, K., Vandemark, D., A uni- 
fied directional spectrum for long and short wind-driven waves, Journal of Geophysical 
Res. vol 102, no. C7, p. 15,781-1,796, July 15, 1997. 

[Garrison et al., 1996] Garrison, J.L., Katzberg, S.J., Howell, C.T., Detection of Ocean 
Reflected GPS signals: theory and experiment, IEEE Southeaston '97, Blacksburg, 
USA, April 12-14, 1997 

[Garrison et al, 1999] Garrison, J.L., Russo, A., Ferebee, M.J., Mickler, D., Armatys, 
M., The GPS Ocean Reflection Experiment on Sapartan 251, ION-99 (Institute of 
Navigation of US), Nashville, Tennessee, Sep 14-17, 1999. 

[Feynman, 1985] QED, the strange thery of light and matter, Princeton U. Press, 1985. 



54 



[Herring, 1990] Geodesy by Radio Interferometry: Tha Application of Kalman Filtering 
to the Analysis of Very Long Baseline Interferometry Data, J. Geophys. Res., Vol 95 
No B8, 12,561-12,581, August 10 1990. 

[Komjathy et al, 1998a] Komjathy, A., Zavorotny, Axelrad, P., Born, G., Garrison, J., 
GPS signal scattering from sea surface: comparison between experimental data and 
theoretical model, presented at the 5th International Conference on Remote Sensing 
for Marine and Coastal Environments, San Diego, 5-7 October, 1998, submitted for 
publication to Remote Sensing of Environment. 

[Martm-Neira, 1993] Martm-Neira, M., A passive reflectometry and interferometry system 
(PARIS): application to ocean altimetry, ESA Journal, vol 17, pp 331-355, 1993. 

[Ruffini et al, 1999] Ruffini. G, E. Cardellach , A. Rius, J.M. Aparicio, Report WP1000 
to ESA Contract 13461/99/NL/GD, Utilization of Scatterometry Using Sources of 
Opportunity. 

[Thomson et al, 1990] Thomson, R.A., Moran, J.M., Swenson, G.W., Interferometry and 
Synthesis in Radio Astronomy, Wiley and Sons, 1990. 

[Zavorotny et al, 1999] Zavorotny, V.U., Voronovich, A.G., Scattering of GPS signals from 
the ocean with wind remote sensing application, IEEE Trans. Geosci. Remote Sensing, 
in press. 



55 



LI pti 



an* 




] 10 JI'J3W 3Q 
11 m», attends 



Magnitude L1 




D 10 ZD H 40 ED 
Hmt r BBcamta 



FTT<flBJd> 





Tim iv locanda 



Magnitude L2 




D 10 10 50 40 ED 



FFT<flBld> 





C la 20 3D « 30 
Dm*, irojndo 



Magnitude L5 




D- 10 £0 3D 40 ED 
Tlmn. -sKDmta 



Fn<fiBiti> 





D 10 Z0 3D « 3D 



Magnitude L12 




1.0 
OJJ 
fa OA 

T 

a 

0.4 
M 



G- 10 20 3D 40 ED 
TTmn. -sscanfta 

ITT(flBld> 




D 10 ZO 3D « 3D 



Hagnitudo L2S 




D 10 20 3D 40 ED 
Titan, -sacond] 



FFTtflekD 




* 1O3JOJ04O50 

Hi 



Hi 



Hi 



103t'30«5C' 

Hi 



1O»30«5O 

Hi 



